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A Mexican Hat with holes: calculating low resolution 
power spectra from data with gaps 
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1 INTRODUCTION 

The calculation of the Power spectrum through direct 
Fourier transform of 2D data in astrophysics is often ham- 
pered by two problems. 



I. Images may have irregular boundaries or parts of the 
image are missing. In many cases masks are applied to re- 
move contaminating foreground sources, creating holes in 
the data that are hard to correct for in Fourier analysis. 
Such problems arise, for example, when angular fluctua- 
tions of a diffuse emission are analyzed and a number of 
compact sources hav e to be excised from the image (e.g. 
IChurazov et~aill2012[ ). 

II. Another problem often encountered when dealing 
with limited data sets is the presence of large scale struc- 
tures, which are not fully covered by the image. The large 
scale power can leak into the observable Fourier frequency 
range, distorting the measured spectrum. This often occurs 
in ID timing analysis when the noise process monitored has 
power on t imescales longer t han the total length of the time 
series (e.g. IScott et al.ll2003h . or in 3D data cubes of hydro- 
dynamical simulations for example, when characterizing the 
turbulent ve locity field in a sub - volume of a larger simulated 
volume fe.g. lDolag et alj|2006h . 



ABSTRACT 

A simple method for calculating a low-resolution power spectrum from data with 
gaps is described. The method is a modification of the A- variance method previously 
described by Stutzki and Ossenkopf. A Mexican Hat filter is used to single out fluctu- 
ations at a given spatial scale and the variance of the convolved image is calculated. 
The gaps in the image, defined by the mask, are corrected for by representing the 
Mexican Hat filter as a difference between two Gaussian filters with slightly differ- 
ent widths, convolving the image and mask with these filters and dividing the results 
before calculating the final filtered image. This method cleanly compensates for data 
gaps even if these have complicated shapes and cover a significant fraction of the data. 
The method was developed to deal with problematic 2D images, where irregular de- 
tector edges and masking of contaminating sources compromise the power spectrum 
estimates, but it can also be straightforwardly applied to ID timing analysis or 3D 
data cubes from numerical simulations. 

These problems are illustrated in the left panel of Fig[TJ 
which shows the Fourier power density spectrum (PDS) cal- 
culated for a 2D image. As input an image with steep PDS 
oc fc -11 / 3 was used. The red points show the Fourier PDS 
for the whole image, blue points show the PDS calculated 
for a section one third of the linear size of the original image 
and green points show the case when about 25% of the data 
in the original image are missing and replaced with zero. 
Changes in the slope and normalization are readily visible 
for blue and green curves. 



A practical way to deal with th ese problems has bee n 
suggested in a serie s of papers by Stutzki et all (|l99ST ) : 
iBensch et all (|200lJ ): lOssenkopf et all (|2008h . Their A- 
variance method preferentially selects fluctuations at a given 
spatial scale a by convolving an image with two filters - a 
compact "core" filter and a more extended "annulus". The 
core filter has a characteristic size ~ a and both filters are 
normalized to unity. The difference between the images con- 
volved with these two filters is an image where all fluctu- 
ations with sizes much larger or much smaller than a are 
suppressed. Therefore the variance of the resulting image 
provides a measure of a typical amplitude of fluctuations of 
size ~ a in the original image. Introduction of a mask helps 
to deal with the boundaries or data gaps. 
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In this paper we discuss a simple method, based on 
the A-variance approach to compute a low resolution power 
spectrum (strictly speaking the convolution of the true 
power spectrum with a broad filter) even when a large frac- 
tion of the original data is missing. The right panel in Fig. [T] 
shows power spectra computed with our method for the 
same images used in the left panel. Compared to the origi- 
nal A-variance method, our approach uses a different imple- 
mentation of the filter, which simplifies the procedure, while 
cleanly compensating for missing data. We further test this 
particular implementation for a number of potential appli- 
cations - e.g. measuring the power spectrum of the surface 
brightness fluctuations of the X-ray images of galaxy clus- 
ters, or characterizing the power spectrum of a turbulent 
velocity field in simulations. 

Many other methods have been devised to deal with 
data with gaps. Notably, multi-taper analysis techniques 
have been develo ped to e s timat e the power spectra of 
one-dimensional ( Thomson! [l982) and multidimensional 
l| Alfred fc Hanssenll 19971) data in Cartesian coordinates and 
on the sphere ( Wieczorek fc Simondl2007h . The choice of ta- 
per functions is adapted to the type of data to be treated 
and the number of tapers is chosen to balance the bias and 
variance in the resulting power spectrum. These methods 
can recover very high resolution power spectra, suppressing 
the power leakage produced by the finite data length and 
are ideally suited to well sampled data where only few or 
no points are missing. Severe gaps in the d ata complicate 
the ap plication of these methods however. iFodor fc Starkl 
()l998h deal with large gaps by calculating separate power 
spectra for each segment and averaging together the results, 
which is not easily applicable to data sets wi th varying data 
and g ap lengths or data in more than ID. IFodor fc Starkl 
(2000) uses multi-tapering methods to compute the power 
spectrum of complete data sets with few small gaps, which 
requires the calculation of optimized taper functions for the 
precise structure of gaps in the data. Although this method 
does a very good job at recovering the power spectral shape, 
it requires long and complex calculations which are not di- 
rectly transferable between data sets with different sampling 
patterns and was only tested for cases where gaps cover a 
small fraction (e.g. 5%) of the data. 

The method we discuss in this paper is simple, ro- 
bust and computationally fast. It has no tuning parameters 
and the interpretation of the resulting power spectrum is 
straightforward. Most importantly, it can be applied with- 
out modifications to the data severely affected by gaps of 
different sizes. The trade-off is its low spectral resolution, 
so it is useful for cases when the underlying power spec- 
trum is a smooth function of a wavenumber/frequency. Pos- 
sible applications in astrophysics include aperiodic variabil- 
ity patterns normally found in AGN light curves; analysis 
of fluctuations in 2D images, e.g. maps of molecular lines, 
X-ray images, Faraday Rotation Measure maps defined in 
irregularly shaped regions; characterization of 3D density or 
velocity fields in numerical simulations. All these cases are 
often affected to a varying degree by gaps in the data. These 
gaps arise mainly due to time constraints in the observations 
in time series, co-adding of several 2D images with different 



orientations and excision of contaminating sources and lim- 
ited computational volumes in 3D simulations. 

We demonstrate that for data sets of sufficiently large 
dynamic range the proposed method recovers well the overall 
shape and normalization of the spectrum even when gaps 
occupy large fraction of the data set. 

In terms of uncertainties in the power spectrum esti- 
mation at a given frequency, the method is analogous to the 
regular Fourier power spectrum, for data sets without gaps, 
provided that similar binning of the Fourier powers is made 
in the frequency space. 

In the following section we will describe the method 
(Sec. [5| and demonstrate its ability to recover the origi- 
nal power spectrum from simulated 2D images with gaps in 
Sec. [3] We also apply the method in spherical coordinates 
and use it to evaluate the power spectrum for a character- 
istic case of CMB data analysis in Sec. 13.21 We explore the 
properties and applications of the method to 3D data cubes 
in Sec. |3] and ID time series in Sec. [5] Expected scatter and 
errors are discussed in Sec. [6] and we summarize our conclu- 
sions in Sec. [7] 



2 METHOD 

We consider an isotropic homogeneous random field, so that 
it can be characterized by a power spectrum P(k), where k 
is a scalar. Instead of calculating the "true" power spectrum 
P(k), where k is the wave number, we want to evaluate the 
amplitude of fluctuations for a broad interval of wave num- 
bers Afc ~ k. We allow for gaps in the data and want to 
recover the correct shape and normalization of the power 
spectrum, provided that P(k) does not contain sharp fea- 
tures. 

The method consists of few simple steps for any given 
spatial scale a: 

(i) the original image is convolved with two filters 
(e.g. Gaussian) having different smoothing lengths o\ = 
cr/s/1 + e and ai = cry/1 + e, where e< 1. 

(ii) convolved images are corrected for the data gaps (see 

32HJ). 

(iii) the difference of two images is calculated. This differ- 
ence image is dominated by fluctuations at scales ~ a. 

(iv) The variance of the resulting image is calculated and 
re-casted into an estimate of the power. 

The variance values are collected as a function of length scale 
or, correspondingly, wave number k r oc l/cr to produce the 
power spectrum P(k r ). The tilde is added to distinguish the 
result from the true power spectrum. 

We first consider the case of data without gaps ( H2.ll 
steps i,iii and iv) and then discuss the procedure of correct- 
ing for gaps ( i)2.2p using a Mexican Hat filter as an example. 



2.1 Data without gaps 

To isolate structures of a characteristic length-scale, we first 
smooth the image / with two Gaussian filters (for simplicity 
we consider ID case): 
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Figure 1. Left: Fourier power density spectrum (PDS) calculated for a 2D image. As input a 351 X 351 pixel image with steep PDS 
(oc fc" 11 / 3 ) was used (shown by the black line). The spectra were obtained by averaging the power spectra, recovered from 50 random 
realizations of images with the same power spectrum. The errorbars were estimated from the scatter between recovered PDS. The red 
points show the Fourier PDS for the whole image, blue points show the PDS calculated for a small section (115 X 115 pixel) of the original 
image and green points show the case when about 25% of the data in the original 351 X 351 pixel image are missing and replaced with 
zero. Clearly both the normalization and the shape of the PDS are modified when parts of the data are missing or significant power is 
present on spatial scales larger than the size of the image (case of image subsection). Right: Power estimated using the Mexican Hat 
filter, corrected for the gaps in the data. The same set of data is used as in the left panel. Clearly the result is insensitive to the presence 
of the data gaps. There is a weak bias in the normalization, which is discussed in the text and in the Appendix FBI 



G a {x) = 



(27T(7 2 ) 1 /2 



(1) 



of slightly different widths <n = cr/vl + e and 02 = 
ay/1 + e, where e < 1. The top panel in Fig. [2] shows an 
example of two such Gaussian filters, normalized to unity, 
together with their difference. For clarity only, in the Fig- 
ure we use e ~ 0.25. After convolving the image with each 
of these filters, both resulting images 7i and I2 will retain 
structures larger than a and lose structures smaller than a 
so the difference image Ji — I2 will predominantly contain 
structures with the characteristic scales ~ a. 
For e — > the resulting filter 



F(x) = G ai (x) — Ga 2 (x) oc 



dG a {x) 



(0-2 

P 

x 2 



(Til OC 



(2) 



which is the familiar Mexican Hat filter. Obviously the shape 
of the filter does not depend on e in the limit of e — > 0. In 
practice we use e = 10 -3 . 

In terms of power spectra, the above procedure is equiv- 
alent to multiplying the original power spectrum of the data 
by the power spectrum of the filter, shown in the bottom 
panel of Fig. [2] This panel shows the Fourier transforms of 



the two Gaussian filters, which are of course also Gaussian, 
and their difference. The square of the difference, is the filter 
power spectrum. As shown in Appendix [X] the peak of the 
filter power is at k r — 0.225/er and its width is ~ 1.155 k r , 
independent of e, as long as e <C 1. Here and below we adopt 
the relation between the spatial scale x and a wavenumber 
k in a form k — 1/x without a factor 2-7T. The resulting filter 
is relatively broad and it cannot be made narrower by us- 
ing smaller values of e. This has the implication that sharp 
features in the power spectrum will be smeared out, so our 
aim is to recover the broad band shape and normalization of 
the power spectrum but not narrow features. The choice of 
the Gaussian filters is twofold. Firstly, they provide a good 
balance in terms of the width in real and frequency space. 
The width in the real space is important when gaps or edges 
are present in the data (see next Section) . A small width im- 
plies that the region affected by a given gap does not extend 
over the whole image. Secondly, the use of Gaussian filter 
is computationally convenient, since n-dimensional convolu- 
tion with a Gaussian in a real space can be easily factorized 
into ID convolutions in each dimension. 

Thus the difference between two convolved images, 
which we denote as I(k r ) is 



I{kr) = h - h = G ai * I - G CT2 * / = F * I, 



(3) 
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Figure 3. Top row: simulated image with power law shaped PDS, each subsequent panel shows the result of filtering the image at a 
different spatial scale k as described in the text, from left to right: fe — 1 /L = 0.6, 0.08 and 0.02, where L is the size of the image. Bottom 
row: original image multiplied by the mask, i.e. / X M . Each subsequent panel shows filtered images calculated from the same simulated 
masked image after application of eq. [6]and multiplication by the original mask, i.e. I c (k r ) X M. No spurious peaks or dips are produced 
around the masked regions or near image edges. These parts of the filtered images where M = 1 are used to calculate the variance for a 
given k r . 




0.2 0.4 0.6 0.8 



k=l/x 

Figure 2. Top panel: 2 Gaussians of a ~ 1 used for filtering (red, 
blue) and their difference (black). Bottom panel: Fourier trans- 
forms of the Gaussian functions shown above (red, blue), their 
difference (black) and the difference squared, amplified for clarity 
(magenta). The magenta curve is the filter effectively applied on 
the power spectrum of the image treated, therefore the variance 
of the filtered image is normalized by the area under this curve. 
The transmitted power peaks at k = 0.225/cr. 



can be used to calculate the variance at scales ~ 1/Av and 
evaluate P(k r ) using a simple relation for the normalization 
(see Appendix lA")l . 



2.2 Data with gaps 

We now consider the image with gaps and introduce a mask 
M such that 

_ J 1 where I is defined , . 

| where / is undefined 

Here "undefined" refers not only to gaps, but also to ar- 
eas outside the image boundaries. Essentially for an n- 
dimensional image we treat the whole n-dimensional space 
outside the image boundaries as a data gap. The image I is 
also set to zero in the gaps and outside image boundaries. 
Thus, one can write I — M x Iq, where 7o is the true image 
without gaps, defined over the whole n-dimensional space. 

Direct application of the filter F described by eq. [2] to 
the image / with gaps will produce many spurious struc- 
tures, which are difficult to correct for. However one can 
use the fact that the I(k r ) image can be represented as the 
difference of two smoothed images. Consider, for instance, 
ii = G C r 1 * I. Convolving the image with gaps with a Gaus- 
sian will still produce spurious features, but their amplitude 
can be drastically reduced by dividing I\ by the mask M, 
convolved with the same Gaussian, to produce a corrected 
image h iC : 

T h Gg 1 * (M x Jp) 

v - Mi - G^*M ■ ^ 

Note, that we make the convolution in the infinite n- 
dimensional space without assumption of the image peri- 
odicity outside the image boundaries. Intuitively the effect 
of the division by the convolved mask is clear: the amplitude 
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of ii = Ga\ *I = G lT1 * (M x 7o) is going to be lower close to 
the gaps or close to image boundaries. The convolved mask 
Mi — Gen * M largely shares these properties and the ratio 
j4- will lack an obvious trend of the amplitude decrease near 

h 

the gaps. The same argument applies to 72, c = 77- • Finally 

M2 



h(k r ) = (h,c-h,c)xM = 



G CT1 */ 



Gq 



G a 



M Gq 



M 



xM.(6) 



The final step of the variance calculation is done only 
for the part of the corrected image I c (k r ) where the mask 
M = 1 (see Appendix IX)) . In summary, the biases intro- 
duced by the shape of the image boundaries and holes in 
the mask are corrected for by subjecting the mask to the 
same smoothing procedure and then dividing the smoothed 
image by the smoothed mask. Therefore, the flux lost by 
smearing the image out of the mask boundaries is compen- 
sated by an equal loss of mask area around the edges. When 
approaching a gap, the convolved mask G ai * M decreases 
smoothly from a value of 1 before the gap to well into 
the gap, provided that the gap is larger than the filter size. 
Therefore the convolved mask typically has a value of or- 
der 0.5 or larger at the edge of the original gajQ. Since only 
unmasked parts of the image are used for the variance calcu- 
lations, points where the denominators in Eq. [6] vanish are 
automatically discarded. 

Our approach is analo go us to the A-variance 



method of IStutzki et alj (|l998l ); iBensch et all I^OPl! ) and 
lOssenkopf et aL I d2008T>. In particular , one of the suggested 



forms of the filters in lOssenkopf et all (|200&f ) (see their equa- 
tion 11) leads to a convolution of an image with a Mexican 
Hat filter, in the limit of their parameter v — > 1. In other 
words the difference between filters is the s ame as used here. 
Howev er the individual filters are different. lOssenkopf et al.l 
(2008) use a Gaussian as a "core" filter and a ring-like shape 
function produced by a linear combination of two Gaussians, 
as a "annulus" filter. In our implementation of the Mexi- 
can Hat filter the role of the core and a nnulus filters (see 
IStutzki et al.in~998i : lOssenkopf et af1l2008l for definitions) is 
played by two Gaussians with slightly different widths. Far 
from image gaps and boundaries only the difference between 
the core and annulus filters matters. But in the presence of 
gaps or boundaries the functional forms of both filter starts 
to be important. Representing the Mexican Hat as a differ- 
ence between two Gaussians simplifies the whole procedure 
when correcting for the missing data and ensures that the 
effect of gaps is almost identical for 7i jC and 72, c- 

Figure [3] shows a simulated image with an isotropic 
power law power spectrum of slope a = —2. The top row 
in this figure shows the image decomposition into compo- 
nents of different spatial scales. The bottom row shows the 
same simulated image to which an arbitrary mask is applied 
(left). Other images in the bottom row show the decompo- 
sition performed on the masked image. Clearly, the method 
is quite insensitive to the presence of the mask, as seen in 
Fig. 13 Indeed, the fluctuations on different spatial scales 



1 Unless we are dealing with an isolated data area, small com- 
pared to the size of the filter and surrounded by gaps. 





Figure 4. Mask 1 (left) and mask 2 (right) applied to simulated 
2D images before computing their P(k r ) spectra. The black areas 
represent the discarded parts of the input images. Data gaps in 
real astronomical images are usually less extreme, compared to 
mask 2. We use this extreme example to test the ability of the 
method to recover the input power spectrum even for images this 
heavily affected by gaps. 



can be recovered in the masked image without introducing 
spurious structures near the edges or gaps. 



3 TEST ON SIMULATED IMAGES 

To test the ability of the method to recover the original 
power spectrum of an image, we generated a set of 2D Gaus- 
sian random fields. We simulate the effect of irregular edge 
shapes and excision of contaminating sources by masking 
out regions with different shapes and size scales. The mask 
in Fig. [3]is typical for wide field images when different expo- 
sures are combined and contaminating sources are excluded, 
creating irregular edges and holes in the mosaiced image. For 
the tests below we use more complicated masks, shown in 
Fig. [4] to verify the performance of the method in more ex- 
treme cases. The tests have also been run using the simple 
mask shown in Fig. O recovering the input power spectra at 
least as accurately. 

We checked for overall biases in the power spectrum 
shape and normalization by generating images with ran- 
domized phases and amplitu des to simulate a Gaussian field 
fe.g. iTimmer fc Konidll995T ). Images were constructed with 
power law power spectra of slopes between and -4 and, 
for each case, 100 realisations were produced and masked 
by the patterns shown in Fig. 2] 

Figure [5] shows the variance of the filtered images as a 
function of their characteristic filter scale k r , for different 
slopes of the power spectrum. Simulated images had a size 
of 1024 2 pixels and a central 300 2 pixels section of each 
image was used to calculate the power spectrum. The P(k r ) 
spectra were computed for each simulated image and their 
average is plotted in Fig.[5]using black dots for the unmasked 
images, green squares for mask 1 and blue stars for mask 2. 

A very small deviation from a perfect power law shape 
is apparent at the high frequency end of the recovered power 
spectra even in the unmasked case. This happens because as 
the filter width approaches the size of the resolution element, 
the filter is under-sampled and fluctuations at those scales 
are less accurately recovered. There is a small decrease of 
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Figure 5. Recovered power spectrum of simulated images. The 
markers show averaged P(k r ) spectrum over 100 realizations of 
simulated images. Black dots show the spectrum computed for 
unmasked images. Green squares and blue stars correspond to 
images masked using patters 1 and 2 (see Fig. [4). The slopes of 
the input power spectra are given in each panel. The input spectra 
are normalized to unity at k = 1. 



power, of not more than 10% of the expected value, at k=0.3 
pix _1 . Below this frequency the recovered power spectrum 
is not affected, so if the range of frequencies covered is large, 
the high frequency part can be discarded and the spectral 
slope can be accurately measured. Similar effects are seen in 
ID and 3D. 

The high frequency points in Fig. [5] have the smallest 
scatter since they are averaged over many Fourier modes. 
Therefore, the statistical uncertainty, discussed in Sec. 6 
is small, and can be smaller than the deviations described 
above. To avoid assigning too much weight to these points 
when fitting the measured power spectra, a systematic error 
of 10% was quadratically added to the estimated statistical 
uncertainty, associated with each point. 

The P(k r ) spectra were measured for each masked and 
unmasked simulated image and the average of each set was 
fit with a power law model. The input and recovered spectral 



input 
slope 



No mask 
slope nori 



Mask 1 
slope noi 



Mask 2 
slope 



norm 
norm bias 



0.0 


-0.02 


1.00 


-0.03 


1.03 


-0.24 


1.33 


1 


-1.0 


-1.00 


0.92 


-0.99 


0.94 


-1.04 


1.04 


0.95 


-2.0 


-2.00 


0.94 


-2.00 


0.95 


-2.01 


0.99 


1 


-3.0 


-3.01 


1.12 


-2.97 


1.19 


-2.99 


1.24 


1.25 


-4.0 


-4.05 


1.91 


-4.05 


1.96 


-4.00 


2.53 


2 



Table 1. Results of power law fits to P(k r ) spectra averaged over 
100 realizations of Gaussian random fields in 2D with and without 
the mask. The input power spectra were normalized to 1 at k = 
1 and have slopes ranging from to -4. For the two groups of 
columns under the headings Mask 1 and Mask 2, the masks shown 
in Fig. [4] were applied to the images before calculating P(k r ). For 
each case, a recovered power law slope and normalization are 
given. Deviations in the normalization are largely explained by 
the normalization bias calculated in Appendix iBl 



slopes and normalizations, obtained by fitting a power law 
to the P(k r ) spectra are given in Table [T] 

The spectral power law slopes are recovered with a de- 
viation not larger than 0.24 for any slope probed and not 
larger than 0.05 for slopes of -1 or steeper. The normal- 
ization deviations from unity follow closely the bias expec- 
tation given in Appendix [B] except for masks similar to 
Mask 2 and very flat slopes. This mask represents a very 
extreme case, containing many small scale gaps that distort 
the power spectrum more than larger gaps as is evident in 
the top panel in Fig. [S] and in Table [T] For simpler masks 
and/or steeper power spectra the method recovers the power 
spectral parameters accurately and the effect of the mask is 
negligible. 

The distortions produced by small gaps covering a large 
fraction of the image (> 50%, similar to Mask 2 in Fig. [5} 
is not easy to predict as they depend on the gap structure. 
For these cases, when the measured power spectrum slope 
is flatter than -1, we recommend estimating the distortions 
through Monte Carlo simulations using known input power 
spectra and the same gap distribution as in the real data. 

The accuracy with which power spectral parameters can 
be recovered depends on the amount of data available, since 
this determines the range of scales covered and the density 
of independent Fourier modes, which affects the errors. As 
an example we fitted the power spectrum of each individ- 
ual realization of the images discussed above to measure 
the scatter in the recovered power law parameters. In these 
fits, the slope and normalization were fitted simultaneously. 
The RMS of the recovered parameters are quoted in Table 
[5] The scatter in the recovered slope and normalization in- 
crease for steeper slopes. The scatter is slightly larger for 
masked images, although there is no significant difference 
between Mask 1 and Mask 2. Notice that these RMS val- 
ues characterize the scatter around the biased mean values, 
similar to those quoted in Table [1] 

3.1 Comparison with the A- variance method 

As mentioned in fjl] and 12.11 the method described 
above is a further development of the A-variance 
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No mask 


Mask 1 


Mask 2 


input slope 


slope 


norm 


slope 


norm 


slope 


norm 


0.0 


0.01 


0.02 


0.02 


0.03 


0.02 


0.03 


-1.0 


0.01 


0.02 


0.02 


0.03 


0.02 


0.03 


-2.0 


0.02 


0.03 


0.03 


0.02 


0.02 


0.02 


-3.0 


0.03 


0.05 


0.04 


0.08 


0.03 


0.06 


-4.0 


0.07 


0.18 


0.07 


0.20 


0.05 


0.20 



Table 2. Root-mean-square scatter of recovered power law pa- 
rameters for sets of 100 simulated images for each spectral slope 
and mask. The images are (300 pix) 2 and Masks 1 and 2 are 
shown in Fig. l4l 



method of IStutzki et all (|l99St ); iBensch et aTJ i|200ll ) and 
lOssenkopf et al.l ( 20081 ) . For the periodic data without gaps 
both methods produce mathematically equivalent resultqj. 
We now proceed with a more detailed comparison of the 

filters performance for non-perio dic data with gap s. 

The functional form used in lOssenkopf et alj (|2008l ). in 
the limit of their parameter v — > 1, is as follows : 



4 
4 



(1/2)" 



tt/ 2 (1/2) 2 



(-) 



where I is the size of the filter. It is obvious that the dif- 
ference between these two filters is the Mexican Hat filter. 
However the shape of the second filter Fi iS , nn (x) is very dif- 
ferent from the first one Fi iCOTe (x). Therefore it is expected 
that the impact of edges and gaps will be different for the 
filtered images corrected for these gaps and edges. This is in 
contrast with the representation of the Mexican Hat filter 
as the difference between two Gaussians with almost equal 
width, which guarantees that the gaps have essentially iden- 
tical impact on both Gaussian convolutions leading to the 
filtered images. 

To verify the above conjecture we made a number of test 
with different power spectra and masks. An example of an 
image and its masked version is shown in Fig. [6] As before 
the non-periodic image was obtained by cutting a section 
of a larger periodic image, in which a random realization of 
a Gaussian random field with slope -2 was generated. The 
spectra recovered from the original and masked images for 
both methods are shown in Fig. [7] The red and green points 
were obtained using the procedure described in i]2.2l for the 
original and masked cases, respectively. Both for the orig- 
inal and for masked images the spectrum agrees well with 
the input spectrum, shown by the horizontal dashed line. 
The spectrum, obtained using Eq. [7] shows larger deviations 
from the true input spectrum, especially for the masked im- 
age, shown in black. Similar results were obtained for other 
types of masks probed. We therefore conclude that for non- 
periodic data sets with gaps the filtering based on two nearly 



2 Provided that during filtering the images are treated as periodic 
data sets 





Figure 6. Original image with power spectrum with slope - 
2 (left) and its masked version (right). This non-periodic image 
was obtained by cutting a small section of a larger image. These 
images arc used to compare the performance of the two filtering 
methods. 
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Figure 7. Comparison of the recovered power spectra using the 
difference between two Gaussian and Eq. [7] The spectra were 
multiplied by fc 2 , so that the input spectrum is flat in this plot 
(shown by the horizontal dashed line). The spectra were aver- 
aged over 10 independent realizations. The red and green points 
were derived from the original and masked images respectively, 
using the filter described in this work. The blue and black curves 
correspond to the original and mas ked images respe c tively , pro- 
cessed using the pair of filters from lOssenkopf et al,l l(2008j). For 
non-periodic data sets and/or for data with gaps the filter based 
on two nearly identical Gaussians performs better. 



identical Gaussians provides more robust and accurate re- 
sults. We note that the value of v recommended by the 
authors is v = 1.5. Repeating the test for this value of v 
gave similar results to the v — 1 case and our approach with 
similar Gaussian still proved more accurate. 
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3.2 The 2D sphere: application to CMB analysis 

The same method can be constructed in spherical coordi- 
nates. As an example, in this section we apply the filter 
formalism to the study of temperature anisotropies of the 
Cosmic Microwave Background (CMB) radiation. Let t(n) 
be the CMB temperature in the direction n. A decomposi- 
tion of the t(n) map on the sphere in a basis of spherical 
harmonics is, 



t(n) = ^2 X] ai,mY^ m (n) 



(8) 



1=2 m=— ! 



where Yj, m (n) is the spherical harmonic of order I, m evalu- 
ated at the direction n and a;, m are multipole coefficients. 
We assume that t(n) is a Gaussian isotropic and homoge- 
neous random field. The angular power spectrum Ci for mul- 
tiple I can be estimated as 



cii, 



(9) 



where the averaging is done over m. Figure [8] shows the 
theoretical prediction for the angular power spectrum of 
the CMB temperature anisotropies for different angular 
frequencies (or multipoles I). The units of the y-axis are 
2?, = i(J + l)C,/(27r). 

Similarly to the Cartesian 2D case described in H2.2I 
one can define a filter as a difference between two Gaussians 
defined on a sphere and allow for the data gaps. This yields 
an image (cf. eg. IA9|) 



t a (n) 



G an * t 



M Go 



M 



x M(n) 



(10) 



dominated by a particular angular scale ~ a. As before the 
variance of the resulting image is calculated and divided 
by the power of the filter. In the multipole space the cor- 
responding filter has the same properties as the 2D filter 
in Cartesian coordinates, described in the Appendix IA1 (see 
eq. I A6f) . Namely, 



-l(l+l)af/2 



-l(l+V)ai/2 



d(l + l)a 2 e- !(i+1)CT /2 , 



(11) 
(12) 



where <j\ = cr/yl + e, <T2 = oVl + e and e< 1. 

The expected shape of the power spectrum can be cal- 
culated as: 



D,C,|F g (Ql 2 (2f + l) 
£,|£,(0l a (2i + l) ' 



(13) 



where l r is defined as the multipole where |-F CT (0I reaches its 
maximum. The corresponding angular power spectrum Ci r 
is shown in Figure [8] with the dashed line. Clearly, the filter- 
ing procedure smears out small scale features, but recovers 
the overall shape and normalization of the power spectrum. 
This smoothed power spectrum is what we want to recover 
from the data using the proposed method. 

In real data, the estimation of the angular power spec- 
trum must deal with the pixelized data and with the pres- 
ence of the Milky Way and other contaminants that make 
necessary the exclusion of some pixels on the sky map. We 



used HEALPbJfl to deal with the pixelized maps. The power 
spectrum is evaluated from the variance of the filtered im- 

(tl) as 



ages 



C, = 



(tl) 



Ei^\Mi)\ 2 \Wi, P ,\ 



(14) 



is the window multipole function for the 
HEALPix pixel. 

The presence of a mask usually biases the estimation 
of the power spectrum multipoles (Ci-s) and, at the same 
time, it couples these otherwise independent quantities. To 
test the method we impose a sky mask that covers ~ 24% 
of the sky, including the Galactic plane and the position of 
bright radio sources. We also consider two different hemi- 
spheres of data separately, defined by an equatorial plane 
perpendicular to the direction (l,b) = (0°, 45°) in Galactic 
coordinates. Therefore, each hemisphere analyzes roughly 
38% of the sky. We apply filters, of width of 0.1°, 0.12°, 
0.22°, 0.35° ,0.42°, 0.5°, 1°, 2°, 3°, 4°, 8°, 10°and 20°. An 
example of a masked CMB realization, filtered at two differ- 
ent scales is shown in Fig. [5] When applying our thirteen 
different filters in the two masked hemispheres, we obtain 
Ci r for each hemisphere, displayed by the filled green and 
blue circles, in Fig. [5] The estimates from the two differ- 
ent hemispheres and for each filter scale a follow closely the 
theoretical prediction provided by the dashed line. 

For a single realization of the sky, the relative scatter of 
power spectrum estimates with respect to the average value 
is known to approximately follow the scaling 



ACi, 
C lr 



(2Z + 1)(AI)/, 



sky 



(15) 



which corresponds to the black solid line in the bottom panel 
of Fig. [5] The symbol / s k y denotes the fraction of the sky 
covered in each hemisphere (around ~36% ), and Al pro- 
vides the effective width, in multipole space, of each filter. 
The blue and green circles show that the estimates of the 
angular power spectrum are scattered around the theoretical 
expectation by an amount that is not far from the theoreti- 
cal prediction. This behavior breaks down at smaller angular 
scales, due to effects related to the finite pixel size, which in 
this case is slightly below 7 arc minutes. 



4 3D DATA CUBES 

When analysing data of numerical simulations, e.g. a veloc- 
ity field in hydrodynamic simulations, one often deals with 
3D data cubes. In this section we extend the method on 
three dimensions and compare our power spectrum results 
with calculations using the conventional Fourier transform. 
As before, we generate isotropic Gaussian random fields with 
a power law power spectrum spectrum, in 3D. In particu- 
lar we consider a red noise process with slope —11/3, cor- 
responding to the 3D Kolmogorov power spectrum, which 
is often assumed when dealing with the turbulent velocity 
field. 

3 HEALPix's URL site: http://healpix.jpl.nasa.gov/ 
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Figure 9. Simulated CMB type maps, masked to remove the Galaxy and foreground point sources, and then filtered at two different 
angular sales. The maps show no evidence for spurious feature near the edges of the mask. 



The power spectrum P(k r ) is calculated analogously to 
the 2D case, by filtering the cube with 3-dimensional Gaus- 
sians. The analysis in Appendix[X]is directly applicable for 
n — 3. As before we assume isotropy in the random field, so 
we are only interested in the power spectrum as a function 
of spatial scale. 

The first question we address is whether the shape of 
the red noise is recovered well when the noise extends below 
the lowest wave numbers set by the size of the cube. To this 
end we calculate the power spectrum for the full simulated 
box and also for smaller "trimmed" boxes, which are cut 
from the original larger cube. The mean of 100 power spec- 
tra evaluated through the variance (see Sj2} for the full box 
and trimmed sections, 1/2 and 1/4 of the original size on a 
side, are shown in Fig. 1101 The la statistical uncertainties 
for all three cases are shown on Fig. [10] with gray shadows. 
The mean of 100 cube realisations and their trimmed sec- 
tions are shown by the red lines on each panel, while the 
scatter is represented by the gray area. Power spectra calcu- 
lated through the conventional Fourier transform are shown 
in Fig. 1111 The input power law power spectrum is shown 
by the dashed black line. Since trimmed cubes are not peri- 
odic anymore the power spectrum recovered through Fourier 
transform is strongly distorted by leakage of power from very 
low frequencies. However, we see good agreement between 
the input and recovered power spectra calculated through 
the variance method. The minor discrepancies are only on 
the smallest and highest wave numbers. 

Both discrepancies are caused by the fact that the value 
of the variance is a convolution of the true PDS with the 
filter in frequency space. As a result at low k the power 
leaks out if the full simulated box is used, the simulated 
cubes effectively have zero power at frequencies lower then 
1/L. This effect goes away if a subsection of the original 
cube is used (see also 2D case in Fig. [TJ. 

We now proceed by considering the impact of data gaps 
on the recovered power spectra. We consider 3D masks cov- 
ering different fractions of the cube volume. Slices of masks 
through the box center are shown in Fig 1121 The first two 
masks are generated randomly and the fraction of missing 
data is 50 per cent in the left panel to 85 per cent in the mid- 



dle panel. Often, when one deals with simulations of galaxy 
clusters, large sub-halos or one of the sub-clusters in merg- 
ing systems are excluded from the analysis. Consequently, 
we generated a third 3D mask that mimics a situation when 
large sub-halo is excluded. The slice of this mask is shown 
on the right panel in Fig. 1121 

Figure [T3] shows the power spectrum of the data with 
gaps calculated through the conventional Fourier transform 
and with the variance method. Clearly, the direct Fourier 
method should not be used for data with gaps. The increase 
of the gaps fraction leads to the leakage of power causing 
the flattening of the spectrum and the decrease of its nor- 
malization. At the same time the power spectra calculated 
through the variance are perfectly recovered, with only mi- 
nor changes on small and large wave numbers even in case 
when 85 per cent of the data are in gaps. 



5 ID: TIMING ANALYSIS 

A possible ID application of the method can be found in 
the study of light curves of variable objects. The method is 
of course not suitable for the search of periodicities, but can 
be useful for studies of aperiodic variability, when the broad 
band shape of the power spectrum is of interest. 

The top panel in Fig. [T3] shows a simulated light 
curve typical for long term X-ray m onitoring of an Ac- 
tive Galactic Nucleu s (AGN) (e.g. lUttlev et all |2002| ; 
iMarkowitz et all I2003T ). The power spectra of these light 
curves is normally modeled as a single or broken power 
law, where the slopes and/or bre a k frequencies are of in - 
terest (e.g. iMcHardv etal] l200d : iPapadakis et al. I 120091 1. 
These light curves often have yearly gaps due to visibil- 
ity constraints and different variability time-scales are cov- 
ered by varying the sampling rate. The window function 
produced by these gaps causes spurious features in the 
power spectrum and different approaches have to be used 
to remove these feat ures from the Fourier transform of 
the light curve (e.g. lUttlev et alj |2002|; IMarkowitz! |2010| ; 
lEmmanoulopoulos et aL 20ld ; Kastendieck et al.ll201ll V 



The case of AGN light curves is essentially different 
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Figure 8. Input angular band power spectrum (solid line), the 
result of filtering such power spectrum with scale changing filters 
(dashed line), and the result of computing the angular power 
spectrum with 13 different filters in two different hemispheres of 
a masked single CMB sky realization (green and blue circles). In 
the bottom panel we show the relative error of the power spectrum 
estimates with respect to the theoretical prediction (dashed line 
in top panel) . The solid line provides the theoretical expectation, 
^2/(21 + 1)/ A,/ f s k y with / sky =fraction of the sky not being 
masked. At I values higher than ~ 800 the filter size reaches the 
pixel scale and the errors detach from the theoretical prediction. 
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Figure 10. Power spectra calculated through the variance for 
boxes of different sizes (same as Fig, lilt . Top: initial size of the 
box; middle and bottom: 1/2 and 1/4 of the original size on a 
side respectively Red: the mean of 100 cube realization, gray: the 
range of statistical uncertainties corresponding to 1 a scatter. 



from images and simulated data cubes since normally the 
light curves are not sampled continuously over a regular grid, 
but consist of snapshot observations on irregularly spaced 
time intervals. Long term AGN light curves, as described 
above, suffer strongly from this effect, since the timescale 
coverage is optimized by combining epochs of very different 
sampling rates as shown in the bottom panel of Fig. Q3] For 
our example, frequencies from 3 x 1CF 9 to 3 x 1CF 5 Hz are 
probed using only 1100 light curve points. If it was evenly 
sampled, this range would require 20.000 points. This un- 
even sampling does not pose a problem to the variance cal- 
culation, since the filter convolution can be performed over 
uneven grids. In this case, the mask is simply the same time 
series, with all flux values replaced by 1. Figure [151 shows 
the filtered light curves on different timescales produced by 
this approach, gaps in the time series and varying sampling 
rates do not prevent a clean filtering of the fluctuations. 

In Figure we have only plotted the useful part of 
each filtered light curve, these sections are selected to reduce 
the effect of aliasing, which is potentially a larger problem 
for light curves than for images. In the case of isotropic 
fluctuations in images, fluctuations on length scales shorter 
than the pixel size are averaged out. In the case of light 
curves, unless they are continuously exposed, variability on 
timescales shorter than the sampling interval adds to that of 
the longer fluctuations. This effect aliases power into lower 
frequency bands and distorts the measured power spectrum. 

Aliasing can be diminished by demanding a mini- 
mum number of exposures within the convolution Gaussian 
around a given lightcurve point for this point to count to- 
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Figure 15. Result of the filtering routine on the light curve shown in Fig. 1141 Gaps in the light curve do not affect the filtering process. 
Short timescales are only probed by higher cadence light curve segments. This is achieved by discarding data points that have less then 
6 neighbors within a time space of 4<r for each value of the filter width a while calculating the convolution. The panel on the left shows 
light curves filtered on timescales of 3300, 1100, 366, 121, 13.5, 1.5 and 0.5 days. The panel in the right shows a zoom to the most 
intensive monitoring section, at the beginning of the light curve, for filter timescales of 13.5, 4.5, 1.5 and 0.5 days, from top to bottom. 
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Figure 11. Fourier power spectrum calculated for a simulated 
box of 64 3 . The power spectrum was obtained by averaging power 
spectra calculated for 100 random realizations with the same slope 
a = —11/3. The input power spectrum is shown with the dashed 
line. Blue, red and green curves show the recovered power spectra 
from the initial box size (64 3 cells), and trimmed sections of 1/2 
and 1/4 of the original size on a side, respectively. 




Figure 12. Slices of masks we applied to the data boxes in 
our calculations. Slices are through the center of the box. Black: 
excluded data, M = 0; white: data used for power spectrum cal- 
culations, M = 1. Left and middle panels: masks are generated 
with random gaps. The fraction of discarded data is 50 and 85 
per cent. Right panel: mask mimics the case of excluded sub-halo. 



wards the variance calculation. In other words, if a point in 
the light curve has no neighbors closer than 4a, for a given 
filter width a, then it cannot provide information on the cor- 
responding frequency k r . Therefore, although the point can 
participate in the convolution of the light curve, it should 
not be counted in the variance for that frequency. The proce- 
dure applied on the simulated light curve counts the number 
of neighbors for each point in the light curve while making 
the convolution and later discards points that have less than 
6 neighbors for a given filter width, producing the filtered 
lightcurves in Figure 1151 Therefore, the shorter timescales 
are only probed by the more intensively sampled sections. 

We tested the recovery of the input power spectra for 
unevenly sampled light curves for different minimum num- 
bers of neighbors within 4<r of each point. Figure [T6l shows an 
example of these tests, for simulated light curves of power 
law slope -1. Each power spectrum is the average of 100 
simulations, to average out the fluctuations between differ- 
ent realizations. In the Figure, the black solid line represents 
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Figure 13. Power spectra calculated for data with gaps (one 
realization). Top: Fourier method, bottom: variance method. 
Dashed curves show initial power spectrum. Blue: 50 per cent 
of data is discarded, see the left panel in Fig. 1121 red: 85 per 
cent of data is excluded, the middle panel in Fig. 1121 green: large 
sub- volume of the box is excluded, see the right panel in Fig. 1121 
Fourier method changes the slope and normalization of the recov- 
ered power spectrum, while the variance method recovers power 
spectrum without distortions on a large range of wavenumbers. 
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Figure 14. Top: Simulated light curve, typical of AGN long term 
X-ray monitoring, with varying sampling frequency and contain- 
ing yearly gaps. Bottom: The sampling rate varies with time, 
here the time difference between consecutive observations is plot- 
ted as a function of time, except for the largest gap, after the 
short, high cadence monitoring. 
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Figure 16. The minimum number of neighbor points required 
to count a point in the variance calculation affects the shape 
of the resulting power spectrum. The power spectra above are 
each an average of 100 realizations, requiring neighbors (solid, 
black), two neighbors (red, dashed), six neighbors (pink , dot- 
dashed) and ten neighbors (blue, dotted). Requiring more than 
10 neighbors does not improve the recovered power spectrum and 
reduces the sampled frequency range as explained in the text. 
Error bars denote the RMS scatter in each set of simulations. 



the power spectrum with no minimum number of neighbors 
required, so that all light curve points count towards all fre- 
quency bins. The peak in the middle is the effect of aliasing 
and the drop at high frequencies is the result of counting 
many isolated points in the variance, since these have all 
the same value they reduce the variance artificially. The red 
dashed line uses at least two neighbors and already corrects 
this last problem, while the aliasing is not resolved. Requir- 
ing 6 neighbors (pink dot-dashed line) already solves a large 
part of both problems, while requiring 10 neighbors (blue 
dotted line) produces a power spectrum quite close to the 
input power law. Using more neighbors does not improve 
the recovered power spectrum but it reduces the probed fre- 
quency range, since at high frequencies the sampling rate 
itself limits the number of possible neighbors. 

Finally, the entire power spectrum of the full range of 
frequencies covered by the light curve in Fig. [I4]can be com- 
puted directly. Figure [17] shows the result of applying the fil- 
ter formalism using a minimum of 6 neighbors to simulated 
light curves with the same sampling pattern and a bending 
power law power spectrum. For this Figure, 100 random re- 
alizations were averaged together. At low frequencies the flat 
input slope and extreme mask produces the slope bias ex- 
pected from the results in Table 1 while the high frequency 
normalization is biased upwards by a factor of 1.3, as ex- 
pected from Appendix [B] The approximate slopes, break 
and normalization of the input power spectrum, shown by 
the solid line, are clearly recovered. 

For light curves similar to the one shown in Fig. Q3] 
the power spectrum is usually estimated through Monte 
Carlo simulations of sample light curves with given power 
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Figure 17. Power spectrum of the simulated light curve with the 
same sampling as in in Fig. 1141 with a flat input power spectrum 
which bends to a power law slope of -2 at a frequency 
Hz. The black symbols correspond to P(k r ) averaged over 100 
realizations of the light curve. Error bars denote the RMS scatter 
in each set of simulations. 



spectrum parameters, attempting to recover only power law 
slopes, breaks and normalizations. 

The simplicity of the variance method is that it can be 
applied directly on the full light curve, regardless of the gap 
and sampling distributions and reproduce the overall shape 
with approximately correct normalization of the power spec- 
trum. This will be particularly useful for large samples of 
AGN light curves from large-area repeated surveys, where 
estimation of aliasing and red noise leak for each particu- 
lar AGN becomes impractical. For comparison, power spec- 
trum fitting of the light curve in Fig. Q3] using PSRESP 
l|Uttlev et al.ll2002h takes several hours whereas P(k r ) and 
its error estimate can be obtained in a few seconds on the 
same computer and both methods recover the spectrum with 
similar accuracy. 



6 ERROR ESTIMATION 

As with other methods, it is possible to predict the uncer- 
tainty of the power estimates by assuming certain properties 
of the underlying variability process. In this case we will as- 
sume that the data being treated correspond to a Gaussian 
field, so that independent points in the power spectrum are 
distributed around their mean values as a \ 2 distribution 
with two degrees of freedom, with standard deviation equal 
to the mean. This means that the uncertainty of the power 
for each individual Fourier mode is equal to the power it- 
self. Therefore, to calculate the expected uncertainty in a 
frequency bin, it is only necessary to combine quadratically 
the power at the frequencies that fall within the bin. 

The effect of the filter in frequency domain is to make 
a weighted average of the power spectrum estimates within 
the width of the filter. Therefore, the expected uncertainty, 
E, of P(k r ) can be calculated as a weighted mean of the 
P(k) values that fall under the frequency filter, where the 
values are combined quadratically. 



Independent Fourier modes in each dimension are sep- 
arated in frequency by 1/Ni where Ni is number of pixels 
in dimension i. This spacing needs to be taken into account 
when calculating the uncertainty in order to include the cor- 
rect number of power estimates under the filter. In the case 
of isotropic fluctuations in n dimensions, where we only cal- 
culate the power as a function of k = VE^iOVW) 2 . dif- 
ferent groups of indices producing the same value of k are 
averaged together and this should also be considered in the 
uncertainty calculation. 

The correct density of Fourier modes is easily obtained 
by summing over all the independent modes, or equivalently 
over half the data points, i.e. over points j, from 1 to Ni/2 
in the first dimension and from 1 to Ni in all the following 
dimensions as 



E(k r ) = 



for 1 dimension and 



E{kr) 



E;ii Er=i( p ( fc )^ 2 ,.( fc )) 2 



(16) 



(17) 



for 2 dimensions. For larger n the formula is similar. 

The filter Fk r (k) is centered on k r in frequency domain 
and is given in equation A6. Clearly for larger data sets 
(larger Ni), the independent modes are spaced more closely 
in frequency so more power spectrum points are averaged 
together under the same filter width and the uncertainty 
decreases accordingly as 1/ y/Hf-xNi. Since the number of 
Fourier modes within the filter width grows linearly with 
increasing frequency in each dimension, the expected error 
spectrum is steeper than the powerlaw power spectrum by 
0.5 in 1-D, by 1 in 2-D and by 1.5 in 3-D. As shown in the 
next section, this prediction is well matched by the data. 

The expected error can be calculated in one pass 
through the data points per frequency kr, so it does not 
imply a significant computational effort. 

The assumption of a Gaussian field implies that each 
independent frequency will have its power distributed as a 
X 2 distribution with two degrees of freedom but the average 
over many such points will lead to a Gaussian distribution. 
For low frequencies the number of Fourier modes within the 
filter is small, so the distribution of powers deviates notice- 
ably from Gaussian. However, since the filter is broad in 
frequency, in ID k r — 4/T already contains enough points 
to produce an approximately symmetric distribution of mea- 
sured power. In 2D and 3D cases the symmetry is achieved 
at even lower frequencies since the number of Fourier modes 
within the filter grows as k n . 

The values of P(k) in the error formula can be estimated 
by interpolating between measurements of P(k r ), given that 
the true underlying form is not known a priori. Using the 
measured power spectrum in the error formula automatically 
takes into account the normalization bias and any other dis- 
tortions, so that the resulting errors always represent the 
statistical uncertainty around the measured values. If the 
spectrum is finally rescaled by the normalization bias, the 
errors should be rescaled accordingly. 
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6.1 Comparing expected error with measured 
scatter 

The error formula predicts the scatter of P(k r ) measure- 
ments around the underlying power spectrum for a Gaussian 
noise process. We tested the accuracy of this approach by 
comparing directly the predicted error to the RMS scatter 
in 100 realizations of a noise process. Figure Q15] shows this 
comparison of measured scatter in markers and predicted 
errors in lines, both relative to the measured power spec- 
trum. Solid lines and nearby markers correspond to ID time 
series without gaps with input power spectra slopes of -1 
and -2 (the power spectra themselves are not shown here). 
There is a very good agreement between the measured RMS 
spread and the estimated uncertainty E in the case of contin- 
uous data sets. The error estimate described in this section 
is equivalent to the expected scatter from Fourier analy- 
sis in the case of continuous data sets. The analysis above 
shows that the scatter around the mean power in binned 
Fourier transforms and in P(k r ) spectra are identical, so 
both methods are equivalent when measuring low resolution 
power spectra of data without gaps. 



6.2 Effect of gaps on the error estimate 

Evidently, gaps in the data will reduce the number of avail- 
able independent points and so the uncertainty in the power 
spectrum measurement should increase. As a first step to 
model this effect we scale E by the number of points actu- 
ally available in the data (N a ) compared to the total number 
of points if there were no gaps N: E' = E x y/ N/N a . 

As a further refinement, it is convenient to count the 
number of available points as a function of k r , since small 
gaps do not affect the power estimates at low wavenumbers. 
This estimate is shown by the dashed lines in Fig. 1181 A 
good match to the measured scatter (markers closest to the 
dashed lines) was obtained when points in gaps were dis- 
counted only if the gap was longer than 10% of the filter 
spatial scale (l/k r ). The errors and scatter measurements 
correspond to the same sets of simulated light curves dis- 
cussed above, this time masked by gaps of random lengths 
so that 85% of the points were discarded. The uncertainty 
for each slope increases and the increase is larger at small 
spatial scales, which are affected by more gaps. In any case, 
the increase in uncertainty for any length scale is between 
E and E x \J N/N a , so that for a small fraction of missing 
points the dependence of the error increase on k r can be 
neglected. 

In the case of light curves with uneven sampling and 
very different sampling rates it is more convenient to scale 
the error by the lengths of useful data stretches to the total 
length. As before, useful light curve stretches can be se- 
lected for each frequency k r by requiring a minimum num- 
ber of neighbors for useful data points. All useful segments 
are added to produce T u (k r ) and the error is scaled as 
E x y/Ttot/T u , where T tot is the total length, used to calcu- 
late the frequency spacing. 

When using these error estimates to fit models to the 
power spectrum it is necessary to use only independent mea- 
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Figure 18. Comparison of predicted relative error (lines) with 
measured RMS scatter (markers), for sets of 100 P(k) of simu- 
lated light curves. The solid lines and corresponding markers are 
calculated for continuous light curves of 2048 points and power 
law slopes of -1 (blue lines and markers) and -2 (red lines and 
markers). The dashed lines and markers show the predicted and 
measured errors for the same light curves after 85% of the data 
points have been masked out producing gaps of random length. 
For this range of slopes at least the analytic prediction of the 
uncertainty in Eq. 1171 works well. 

surements that, given the mixing effect of the filter, should 
be separated by at least a factor or 2 in frequency. 

The analytic error estimate was compared to Monte 
Carlo simulations of Gaussian fields in 1, 2 and 3D, for sim- 
ple power law or broken power law power spectra of zero 
or negative slopes. Although this is a limited subset of all 
possible power spectra it does cover a wide range of astro- 
physical phenomena of interest. More examples are given in 
Fig. 1191 Here 2D images with broken power law power spec- 
tra and affected heavily gaps as shown in mask 2 of Fig.Uare 
examined. The decrease in number of points was estimated 
from gap size distribution in the x direction only to estimate 
the errors. Since the gap structure is largely isotropic, this 
simple ID estimate gave a good approximation, shown by 
the dashed line, error estimates from RMS scatter of 100 
trials is shown in black markers. In the same figure we plot 
the predicted error and RMS scatter of 100 simulated data 
cubes with power law slope -11/3, also masking out 50% 
of the data points. The mask in this case is similar to that 
shown in Fig. 1121 which affects all scales in approximately 
the same way so the effect of the mask can be estimated 
simply by rescaling the predicted error by the square root 
of the ratio of the total number of points to the number of 
points actually available. 

6.3 Data with periodic gaps 

A special case that we considered is largely found in 1-D light 
curves where observational constraints produce strictly pe- 
riodic gaps. If the data and gap stretches are of equal length 
or the data stretch is longer, then the variance method and 
error estimate can be directly applied. However, if the gaps 
(T g ) are longer than the data stretches (Td) and the sam- 
pling pattern is strictly periodic, then some timescales are 
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Figure 19. Relative error prediction for 2D (red) and 3D (blue) 
data. In this example the input power spectrum of the 2D data 
is a broken power law of low frequency slope -0.5 and high fre- 
quency slope -2.5, and for the 3D data is a power law of slope 
—11/3. The images have been masked out by random sized gaps 
covering 50% of the image. The predicted errors divided by the 
average power spectrum are shown by the dashed lines and the 
measured RMS scatter of the 100 realizations is shown by the blue 
and red markers. The error formula in Eg. 1171 predicts the errors 
accurately for simple and broken power law spectra in images and 
data cubes as well. 

effectively not probed by the data and the method cannot 
be expected to produce reliable power spectral estimates for 
the corresponding frequencies. Frequencies above (2xTj) _1 
and below (4 x (Td+T B ))~ are properly reproduced and the 
error estimate applies, as shown with an example in Fig. [20] 
. The power spectrum of frequencies within this range are 
not reliable and should be discarded. 

Unsurprisingly, it is possible to recover the power spec- 
trum for the timescales that are well sampled, either shorter 
than each individual data stretch Td or longer than a few 
times the sampling interval. In cases like this it is com- 
mon to make separate estimates for the high and low fre- 
quency power spectrum, by averaging together individual 
power spectra from short stretches for the first case and by 
computing the power spectrum of low resolution light curves, 
averaging together data over Td and using the sampling in- 
terval as time step for the latter. One of the advantages of 
the variance method is that it can compute reliably the same 
range of frequencies directly from the complete light curve, 
without any further manipulation of the data. 



7 CONCLUSION 

A simple method for estimating a low resolution power spec- 
trum from data with gaps is described. Essentially the vari- 
ance associated with a given scale is calculated by convolv- 
ing the image with a Mexican Hat filter. Gaps in the data, 
described by the mask M are corrected for by representing 
the Mexican Hat filter as a difference between two Gaussian 
filters, convolving the image and mask with these filters and 
dividing results before calculation of the final Mexican Hat- 
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Figure 20. Effect of periodic gaps in 1-D data. The above P(k r ) 
spectrum is the average of 100 realizations of slope -2 light curves 
masked by periodic gaps, where every 20 data points are followed 
by 80 gap points. The RMS spread of the simulations is shown 
by the black markers and the uncertainty predicted by Eg. 1171 is 
shown by the corresponding black line. The average power spec- 
trum in blue crosses fits well the input power spectrum shown by 
the top black line and the error estimate works well for frequencies 
above (2 X T d ) _1 = 0.025 and below (4 X (T d + T g ))~ 1 = 0.0025, 
the power spectrum for frequencies in between these values should 
be discarded. 

filtered image. The variance of the filtered image is then 
calculated and the power spectrum is evaluated by repeat- 
ing the procedure for different filter scales. The calculated 
power spectrum is smeared out by the width of the filter, 
so sharp features are lost, but the broadband spectral shape 
and normalization are recovered well. 

The strength of the method is that it is simple and 
robust. It can deal with severe gaps in the data and pro- 
duce accurate power-spectral power law slopes and normal- 
izations with no tuning parameters and is computationally 
cheap. By dividing the filtered images by their correspond- 
ingly filtered masks and maintaining the calculation in the 
space domain, the method cleanly compensates for data gaps 
even if these have complicated shapes and cover significant 
part of the data set. We use simulations to show that the 
power spectrum recovered from complete data sets and from 
their masked versions are consistent and no additional bi- 
ases are introduced by the masks. The method can be ap- 
plied straightforwardly to ID timing analysis, 2D imaging of, 
for example, brightness fluctuations in galaxy clusters, and 
higher dimensional data cubes from numerical simulations. 
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APPENDIX A: CARTESIAN COORDINATES 

Consider a n-dimensional 'image' I(x), where x is n- 
dimensional vector, and corresponding mask M(x) with val- 
ues of either or 1. M(x) = means that this particular 
region of the image does not contain useful information and 
should be ignored when calculating the power. We first con- 
sider Cartesian coordinates. Our goal is to estimate typical 
amplitude of the power spectrum for a given spatial scale a 
or, equivalently, given wave number k r . Here and below we 
adopt the relation k r — l/a without a factor 2n. The image 
is assumed to be an isotropic Gaussian random field so we 
will only compute the power as a function of the scalar k r - 
We start by defining a Gaussian filter in spatial domain, 
which will be convolved with the image: 



G a (x) 



1 



(Al) 



(27TCT 2 )™/ 2 " ' 

Assuming first that the mask is equal to 1 everywhere we 
can combine two Gaussians to create a filter which selects 
fluctuations with a given scale (see Fig. [2]): 



F kr {x) = Gai{x)-G a i{x), 



(A2) 



where o\ = a/\/l + e, <72 = <J\/l + e and e << 1. The shape 
of the filter is identical to the Mexican Hat in the limit of 
small e. In Fourier space the corresponding filter is: 

~ir 2 fc 2 cr 2 _ e -2ir 2 fe 2 o-| (A3) 



F kT (k) = e — " - 1 - e - " " 2 . 
Making a Taylor expansion for e we obtain 
F kr (k)^e4n 2 k 2 a 2 e- 2 * 2k2 ° 2 . 



(A4) 



The filter shape is therefore independent of the value of e 
in the limit of small values so that the Taylor expansion is 
valid. The value of e only determines the normalization of 
the filter and this effect is canceled in the conversion of the 
variance of the filtered image, V kr to P(k r ) below. The peak 

of the filter is at k — — (see Fig. (21). ft is therefore 

natural to relate the characteristic scale k 



to this particular filter width a by: 
1 1 0.225079 



corresponding 



(A5) 



Thus the final expression for the filter in Fourier domain is: 



(A6) 



Convolving the image I(x) with the filter Fk r (x) described 
by ea. lA2l and integrating the square of the convolved image 
is equivalent to integrating the product of the power spec- 
trum P(k) and the square of the Fourier transform of the 
filter, F kr (k). Below, V kr is the variance of the filtered image 
(F * /) since the mean of this filtered image is zero. 

V kr = [ (F * lfd n x = 



P(k)\F kr (k)\ 2 d n k. 



(A7) 



Notice that in eq, IA6l the filter drops to zero for k larger and 
smaller than k r so that the function P(k) can be approxi- 
mated by P(k r ) and taken out of the integral: 



V kr « \ 2 P{k r 

= <L 2 P(k r 



*V e -»(*) ! 



(k r ) 



n 

2 +1 



d n k = 



" 1 7T 2 kr 



(A8) 



The last approximation breaks down for very steep power 
spectra, of positive or negative slope, where its product with 
the filter does not drop rapidly to zero for values of k differ- 
ent to k r . 

The quantity V kr is calculated from the image by follow- 
ing the filter and difference method described and expression 
I A8I relates this quantity to the power P(k) through a scale- 
dependent normalization. Therefore, the power P(k r ) can 
be evaluated. 

In practice the mask is not unity everywhere. In this 
case the above calculation is done by convolving the image 
and mask with the two Gaussian filters, dividing the con- 
volved images; subtracting results and applying the original 
mask to the result: 



S k 



Ga\ *M ~ G~„ 



M 



M 



(A9) 



The square of S kr (x) is then integrated. The regions where 
the image is not defined (S kr (x) = 0), which do not con- 
tribute to the variance, are simply compensated for by 
counting all pixels where M(x) = and making appropriate 
scaling: 



Vfc r ,o6s — 



N 



x / S 2 kr (x)d" 



(A10) 



where N = J d n x and A r (M =i) = / M(x)d"x. 

Comparing eg. IA10l and lA8l we get the final estimate of 
the power density spectrum P kr for a given wave number 

Vk r ,obs _ V kry obs 



P(k r ) = 

where 
T 



: 2 n(f + l)2-t- 1 7rtfc,'J e 2 T(n)fc? ' 



(n) = n(| + l) 2-%-\%. 



(All) 



(A12) 



; , , „ lSjr 3 ^ 2 for „ i 2 , 3. 



and T(n) - 3y ^, .',~ i7r 

The evaluation of the power density spectrum thus re- 
duces to eg. lAlll for a set of values of k r . 



APPENDIX B: NORMALIZATION BIAS FOR A 
POWER LAW SPECTRUM 

For a pure power law spectrum P(k) oc k~ a the expression 
IA7l can be easily evaluated and exact relation between P(k r ) 
and P(k r ) can be written. The shape of the spectrum is 
of course recovered correctly, while the normalization may 
differ slightly, caused by moving P(k) outside the integral 
in eq. IA8I 
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Figure Bl. Bias in the normalization of the recovered spec- 
trum for a pure power law power spectrum, as a function of slope 
for different dimensions of the problem (red - ID, blue - 2D, black 
- 3D) 
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Corresponding bias £ is shown in Fig. IB1I One can see that 
for the most relevant problems (2D or 3D geometry, 4yHr 
in the range [-3:0]) the bias is modest. 
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